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Matrix Product State (MPS) wavefunctions have many applications in quantum information 
and condensed matter physics. One application is to represent states in the thermodynamic 
limit directly, using a small set of position independent matrices. For this infinite MPS ansatz 
to be useful it is necessary to be able to calculate expectation values, and we show here 
that a large class of observables, including operators transforming under lattice translations 
as eigenstates of arbitrary momentum k, can be represented in the Schur form of a lower 
or upper triangular matrix and we present an algorithm for evaluating such expectation 
values in the asymptotic limit. The sum or the product of two such Schur operators is also a 
Schur operator, and is easily constructed to give a simple method of constructing arbitrary 
polynomial combinations of operators. Some simple examples are the variance ((H — E) 2 ) 
of an infinite MPS, which gives a simple method of evaluating the accuracy of a numerical 
approximation to a eigenstate, or a vertex operator (cj^cj^cj^c^). This approach is a step 
towards improved algorithms for the calculation of dynamical properties and excited states. 
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1. Introduction 



The Matrix Product State (MPS) ansatz 0, H forms the basis of many numerical 
algorithms, notably the Density Matrix Renormalization Group (DMRG) 0, 0], 
and Time Evolving Block Decimation (TEBD) [5j. These algorithms can be ap- 
plied directly in the thermodynamic limit of a translationally invariant system 
(invariant under translations of some fixed number of lattice sites) 0-0], which 
gives advantages over traditional finite-size scaling calculations. In the infinite size 
variant of DMRG, the converged fixed point produces a translationally invariant 
MPS, as studied for example by Ostlund and Rommer [9|, which gives a compact 
representation of an infinite size wavefunction, from which the behavior of possible 
correlation functions can be obtained from the spectrum of the transfer operator. 
Imaginary time evolution via the iTEBD algorithm 0] produces the same fixed 
point, and indeed the iTEBD and iDMRG algorithms are very similar [8], the 
main difference being the algorithm for the local update of the tensors in the MPS; 
whereas TEBD uses a local evolution of a single bond, the iDMRG algorithm uses 
a very efficient local optimization of the total energy. A drawback of the iTEBD 
approach is the use of the Trotter-Suzuki decomposition (lo| . which gives a con- 
straint on the size of the unit cell which must be a multiple of 2, and interactions 
beyond nearest neighbor are more difficult, with typical implementations using rel- 



atively inefficient swap gates U|. Even for nearest-neighbor interactions only, the 



computation time is linear in the size of the unit cell. On the other hand, iDMRG 
allows any size unit cell that is compatible with the periodicity of the wavefunc- 
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tion, and the performance of the algorithm is essentially independent of the size 
of the unit cell, and depends principally on the number of local optimization steps 
which can be interpreted as the total size of the lattice (which, for good conver- 
gence, will scale with the longest correlation length in the system). Nevertheless, 
(i)TEBD does have some advantages. As described in 0, Q3] the 2-site DMRG 
algorithm introduces a small perturbation at the location of the bare sites, which 
is a consequence of a non-zero truncation error. This causes the variational state 
to be less than optimal for a given basis size. The effect is small, nevertheless for 
approaches based on scaling with respect to the basis size it is an unwanted effect 
and the remedy is to approach the converged fixed point in a gradual way, which 
corresponds to imaginary time evolution with a small time-step (this can of course 
be performed in DMRG, without a Trotter-Suzuki decomposition, by replacing the 
eigensolver with a multiplication by 1 — eH). For real-time evolution, the time- 
dependent DMRG algorithm [l3l . 14 ] is indistinguishable from (i)TEBD, save for 



an unimportant change in notation. 

For a homogeneous system, an infinite MPS offers many advantages over a finite- 
size MPS. The absence of boundaries avoids many of the problems with Friedel 
oscillations that complicate the calculation of correlation functions in finite-size 
DMRG On the other hand, conventional finite-size scaling with respect to 
the lattice size is not possible, but instead this can be replaced by a scaling with 
respect to the largest correlation length in the system £, which for a critical state 
scales with the number of states kept m in the MPS auxiliary basis with a power 
law fui\v&, 



f = m K . (1) 



The exponent k is a function only of the central charge of the conformal field theory 
describing the critical point [l8( ]. 

The evaluation of local or finite-range expectation values on an infinite MPS is a 
straightforward calculation. On the other hand, in [8j], the general approach for cal- 
culating infinite sums of local terms on an infinite MPS was presented, for example 
to compute the energy per site and the fixed point matrices of the Hamiltonian op- 
erator. In this paper, we extend these results to present, in detail, an algorithm for 
constructing the expectation values and auxiliary matrices of an arbitrary polyno- 
mial function of such operators, which includes operators at non-zero momentum, 
fermionic operators, and string operators. This generalizes the results presented in 
0] for two- and three-body operators. After discussing some background on ma- 
trix product states in the infinite size limit, we describe in detail the algorithm 
for obtaining expectation values of triangular MPO's in Sec. [3l and as a simple 
example of these techniques we discuss the utility of using the variance a 2 of the 
energy as a convergence measure in numerical MPS algorithms in Sec. HI Finally, 
we summarize the results and give some concluding remarks. 



2. Infinite size MPS 

A position independent MPS on an infinite lattice is represented by the form 

J2'--A Sl A S2 --- | | s 2 > <8» ■ ■ ■ • (2) 

{Si} 

The local index represents an element of the d-dimensional local Hilbert space 
at site i € Z of the infinite lattice, and the matrices A s have dimension m x m. 
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Note that sometimes in the literature x or D is used instead of m. In general, we 
need not require that the unit cell of the state is exactly one lattice site, but we 
can have in principle any finite periodicity. Although we present all results here 
for the case of a 1-site unit cell, all of the results presented in this paper generalize 
straightforwardly to the case of a multi-site unit cell, at a computational cost that 
is linear in the size of the unit cell. MPS wavefunctions of this form Eq. ([2]) have 
been studied extensively in the literature, eg in and 0. For an introduction 
to the algorithms for computingthe elements of the infinite MPS representation 
numerically, see Ref. [a, [8(. In [9(, the MPS ansatz was generalized to a simple 
representation for an excited state (Bloch state), which is a generalization of the 
single- mode approximation (SMA), constructed by inserting an additional matrix 
at all possible positions in the lattice, which we write here in the limit of an infinite 
size lattice, 



Q,k)=Y^ eikiA ''''' Aai ~^ A ' iA ' i+ '''' \*i)---\s j - l )\s j )\s j+l )-- - . (3) 



This can be represented in a more compact form as a triangular MPS, closely re- 
lated to the W^-state as a position independent MPS with matrices A s of dimension 
2m x 2m, given by, 



In this formulation the final wavefunction is accumulated in the bottom-left ma- 
trix entry, rather than on the the diagonal. This is a consequence of the explicit 
breaking of U(l) charge symmetry to generate a particle-like excitation, thus a 
finite-dimensional MPS representation that transforms as a U(l) invariant (scalar) 
MPS is not possible [lij . This form is closely related to the application of a mo- 
mentum k triangular operator to a state, for example the matrix product operator 



which constructs a boson of momentum k, where is the boson creation operator 
and / is the identity operator for a single site of the lattice. 

A key advantage of the MPO formulation [20(] compared with the ad hoc methods 
traditionally employed in matrix product numerical approaches, is that the MPO 
easily allows arithmetic operations, making the construction of complex operators 
rather simple. In particular, sums and products of MPO's are constructed simply by 
taking the matrix direct sum and matrix direct product respectively of the MPO 
matrices, and often the resulting MPO can be factorized to reduce the matrix 
dimension. For example, an MPO representation of the number operator Nf~ = 




(4) 



(MPO) 




(5) 




I 1 \ 

&t e ik I 

b e~ ik I 
\6+6 e ik b e~ ik tf I ) 



(6) 



The MPO formulation also allows for a convenient and efficient representation of 
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longer range interactions, as a sum of terms that decay exponentially with distance 



211 ] . A surprisingly small number of such terms can be used to approximate a long- 
range polynomial decaying interaction. For a finite system, the expectation value 
of a triangular MPO operator will be a function system length that may have 
a complicated short-range behaviour. In this paper we calculate the asymptotic 
functional form of the expectation value in the limit as we approach an infinitely 
large system, which approaches a polynomial function. To see how this works, 
consider a simple example of the Hamiltonian operator, for example for the Ising 
model in a transverse field, 

= $>K +1 + A5>f , (7) 



which has the MPO representation 



H= ( a z ) . (8) 

On an infinite lattice this expectation value diverges, with the physically relevant 
quantity being the energy per lattice site. Equivalently, we can consider an expec- 
tation value on n sites of the lattice, which we denote (H) n , which in this example 
is equal to nEo, where Eq is the energy per site. More generally, the expectation 
value of an arbitrary triangular (Schur) operator will be some polynomial in n. 

On a section of the lattice of length n, the exact expectation value of the MPO 
is not well defined unless a boundary condition is specified. For example, with 
open boundary conditions and the Ising Hamiltonian, the erf term occurs exactly 
n times, whereas the nearest-neighbor term erf erf , j occurs only n — 1 times. This is 
a consequence of the boundary condition, and in the asymptotic limit the resulting 
constant term is not relevant and in fact is problematic when constructing higher 
order operators, where such boundary contributions can potentially affect all sub- 
leading terms in the polynomial form. These spurious terms are easy to eliminate 
however, by constructing the correct fixed point equations for the expectation value. 
This procedure is described below. 



3. Fixed point equations of triangular MPO's 

For each index of the M x M dimensional MPO, we can associate a matrix E (in 
DMRG notation, this is called the block operator). To find the expectation value 
of these operators, we can use a recursive formula based on the notion that the 
operators E^ are a function only of the previously calculated Ej, for j > i, because 
of the specific triangular form of the MPO. This gives a recursive algorithm for the 
expectation value, whereby we calculate Em, Em-i, ■ • ., E\ in turn. Solving each 
Ei matrix will require 0(dm 3 ) operations, giving a total computational complexity 
of 0(m 3 dM 2 ). 

We define the transfer operator, which acts on the m x m E matrices, 

T(E) = Y J A S ^EA S . (9) 

s 

This operator has a spectral radius of 1, and the eigenspectrum determines the 
scaling form of all possible correlation functions [jj. If the wavefunction is both 
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parity and time inversion invariant, then we can choose a normalization such that 
the A s matrices are symmetric, in which case the transfer operator is also sym- 
metric. If the wavefunction is invariant under only the combination of parity and 
time (PT), then the A s can be chosen to be Hermitian, in which case the transfer 
operator is Hermitian. For a wavefunction that is only CP or CPT symmetric, 
the A s matrices can be chosen to be symmetric or Hermitian in combination with 
a local basis transformation (corresponding to charge inversion), in which case 
the transfer operator is not Hermitian, but is a normal operator. We assume that 
there is one eigenvalue of T equal to 1. This follows from the normalization of the 



wavefunction [9j, |20J, and corresponds to the left/right eigenpair of the identity 
operator and reduced density matrix. In general there may exist more than one 
eigenvalue equal to 1, which signals long range correlations in the state. For the ex- 
position of the procedure for calculating expectation values, we assume that there 
is only a single eigenvalue 1, but the generalization to long range correlated states 
is straightforward and is discussed in Sec. 13.31 

Generalizing this transfer operator, we can let X be an operator acting on the 
local Hilbert space of the MPS, and define 

T X {E) = J2(s'\X\s)A s ^EA s . (10) 

s's 

Now given an MxM dimensional MPO W, we can express the action of adding one 
site to the expectation value in terms of the polynomial form for the M different 
matrices Ei, 1 < % < M, as 

Ei(n + 1) = T Wii {Ei{n)) +$^2V J4 (E,-(n)) . (11) 

j>i 

The reason why we have split this into two terms, with the diagonal part Wu 
and the off-diagonal part Wij, is that the Ej matrices for j > i are assumed to 
be already calculated, so the off-diagonal part is some matrix function of n. Let 
Ci(n) = Ylj>i (Pj( n )) be the fixed right hand side, and let X = Wu be the 
diagonal element of the MPO, that acts on the local Hilbert space. Then Eq. (jlip 
reduces to 

Ei{n + l) = T x {E i {n)) + C i (n) . (12) 

The operator representing the observable is Ei(n), which has an expectation value 
of Tr pEi(n), where p is the reduced density matrix (the right eigenvector of the 
transfer operator T with eigenvalue 1). To solve these equations in the large n 
asymptotic limit we consider several cases, firstly zero momentum and later we 
generalize to non-zero momentum and string operators. 



3.1. Zero-momentum 

An MPO containing only zero momentum components is characterized by the 
diagonal components Wu being proportional to the identity operator, Wu = xl, 
with the prefactor x satisfying either x = 1 or \x\ < 1. In this case, the C and E 
matrices are polynomial functions of n, with matrix-valued coefficients. Therefore, 
let 
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where C m and E m are matrix-valued coefficients of the p and p + 1 degree poly- 
nomials C(n) and E(n) respectively. For clarity of notation we have suppressed 
the subscript i from the C(n) and E(n) matrices, with the understanding that the 
same form of fixed point equations will be solved for each column of the MPO. Now 
let diagonal operator Wu = xl, where x is a c-number. We can further divide into 
sub-cases: If x = 0, then Eq. (fl~2j) reduces to simply E(n + 1) = C(n), or equating 
coefficients, 

E m = C m - (M) 

j=m+l ^ ' 

which again is obtained recursively starting from E p , . . ., E\. 

The second case is \x\ < 1. In this case, Eq. (fT2j) reduces to 



Ei(n + 1) =xT{Ei(n)) + C(n) . (15) 
Equating coefficients of the polynomial expansion, we get 



[1 - xT)(E m ) = C m - \L) E i- ( 16 ) 

j=m+l 



For each index m = p,p— 1, • ■ • , 0, the right hand side is a fixed matrix so the E m is 
obtained as the solution of a set of linear equations. The solution to these equations 
corresponds to taking the limit n — > oo, such that the geometric series defined by 
Eq. ()12p converges to a fixed point. Since we have \x\ < 1, the operator 1—xT is non- 
singular and the solution is unique. If the transfer operator T is Hermitian, then this 
set of linear equations can be solved using a simple conjugate gradient method. For 
more general cases, conjugate gradient is not suitable but the GMRES algorithm 
gives good convergence at the cost of higher memory requirements, although this 
is typically not a significant limitation in this context. 

The final case for zero- momentum operators is when the prefactor x = 1. Because 
the transfer operator T has an eigenvalue 1, this means that the left hand side of 
Eq. (|16p is singular and it is necessary to decompose the matrices E(n) and C(n) 
into components parallel and perpendicular to the identity (ie. the eigenvector 
corresponding to the eigenvalue 1 of T). We do this via 

Em — &mJ + Em (~\7\ 

C — o T 4- C \ > 

^m — ^m 1 ~r 

where E m denotes the component of E m perpendicular to the identity operator and 
e m is the coefficient of the component of E m in the direction of the identity operator. 
Now, the components E m are well defined by Eq. f)16[) . since these matrices are 
orthogonal to the singular component of 1 — T. In a numerical implementation it 
is necessary to take special care however to ensure that the solution to the linear 
equation is stable, by removing any spurious components in the direction of the 
identity during the course of the solver algorithm. 
The components in the direction of the identity satisfy 



1 

e m +i 



m + 1 



P+i 

E 

k=m+2 



<■„< ^ . )< I, 



(18) 
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which again can be solved straightforwardly starting from m = p down to m = 0. 
Hence if c p ^ 0, then the degree of the polynomial will be increased by 1, as we 
will end up with a non-zero component e p +\. Note also that cq is not defined in 
this procedure, which is what we expect since the constant offset of the expectation 
value is a boundary term. We are free to choose eo = 0, which is the choice that 
removes spurious sub-leading boundary contributions. 



3.2. Finite momentum 

An MPO involving operators of finite momenta introduces oscillating components 
into the expectation value, that cannot be represented by a polynomial. However, 
as long as there is only a finite number of momenta, which will always be true 
for a finite-dimensional MPO, we can write each coefficient matrix E(n) and C(n) 
uniquely as a sum over momenta, as 

C(n) = Y, k e lkn C^(n) 

E(n) = J2 k e ikn E( k \n) { ' 

where the k summation is over all distinct momenta with non-zero contributions. 
We can then expand the C^ k \n) and E^ k \n) as polynomials in n as before, and 
equate coefficients of e lkn n m for each k, m. We now consider each of the possible 
cases for the diagonal operator X = Wa, generalizing the results from the previous 
section. 

Firstly, if X = 0, then Eq. (fT^) acquires a phase factor, giving 

i#> = e -*CjW- £ P( j )^ k) , (20) 

j=m+l v 7 

and similarly, when X = xl with \x\ < 1, we get a phase factor difference from Eq. 
flZD), giving 



(1 



-ik 



xT){E r , 



p~ ik r 



p+1 

E 

j=m+l 



J 



(21) 



Now the third case, for X = xl with \x\ = 1, let x = e . We have two distinct 

(k) 

possibilities. If there is a non-zero component Cm with the same momentum, then 
we get a diverging component and we must again treat the components parallel 
and perpendicular to the identity (eigenvector of the transfer operator T with 

~ (k) 

eigenvalue 1) separately. For the perpendicular components E m , we get 

p+1 , 

(1-T)(^) = e-^)- £ ( 3 WK (22) 

j=m+l ^ ' 
(k) 

while the component in the direction of the identity, e m , is 

(fc) _ 1 
e ' m+1 m + 1 



-ik c (k) 



P+i 

E 

j=m+2 



(23) 
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(k') 

For components C m with a different momentum k f= k, the components perpen- 
dicular to the identity converge, 



(l-eW-VT)(EW) = e- ik 'cW- J2 (!r)^ k '\ (24) 

j=m+l ^ ' 



while the component in the direction of the identity acquires a new component at 
the momentum k, 



(fc/) = C m 



ik') 

(25) 



e ™ ) ~~ e ik> _ e ik ( 26 ) 



(k0 



This applies in an additive sense, so that if there are many different components 
Cm ^ with different momenta k', the coefficient is the sum of all ik f™ ik ■ 

For the finite momentum case, the final matrix E\(n) will contain in general 
many oscillating components with different momenta k. The actual expectation 
value is obtained when every term contributes an integer number of wavelengths 
A = 2ir/k. Therefore, we can restrict n to be the lowest common multiple of all 
wavelengths A, in which case all of the oscillating components e lkn become equal 
to 1 and we recover a simple polynomial. In fact, this procedure also works if the 
A are irrational, in the sense that for a large enough n we can get all components 
arbitrarily close to containing an integer number of wavelengths. 



3.3. String operators 

Now that we have the machinery to handle non-zero momenta, the generalization 
to operators on the diagonal Wa = X which are not proportional to the identity 
operator is straightforward, and in this section we sketch the solution. Applications 
where non-trivial operators occur on the diagonal of a triangular MPO include 
string order parameters, where the diagonal component is some unitary operator, 
and fermionic operators, where the diagonal component will be the local number 
parity operator (— 1)^, corresponding to a Jordan- Wigner string. We can assume 
that it is normalization of the operator is such that the spectral radius of the 
transfer operator Tx is at most 1 (otherwise the expectation value may diverge 
exponentially). We can easily relax the requirement that there is at most one 
eigenvalue of Tx equal to 1, which covers the case where the wavefunction contains 
long range correlations. To solve the fixed point equations for such an operator, we 
must identify the eigenvector subspaces of Tx with eigenvalue 1 and eigenvalues 
of the form e lk , with norm 1. The components of the Ei(n) operators orthogonal 
to these subspaces converge to a fixed point in the same way as we have treated 
previously, according to Eq. (|20l) . Eq. (|21j) and Eq. (|24j) . For each component in the 
direction of an eigenvector of Tx with eigenvalue of norm 1, we treat the coefficient 
in the same way as for the components in the direction of the identity operator, 
Eq. (|23p and Eq. (|26p which must be done separately for each such component. 
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Figure 1. Energy versus variance in the spin 1/2 Heisenberg model, with a number of states kept varying 
from m = 40 to m = 200. Inset shows a closeup of the points used to obtain the linear fit. 



4. Variance as a test of convergence 

A simple application of this method is to calculate the expectation value of the 
square of the Hamiltonian operator. This gives a polynomial of degree 2, and for 
the typical case where there are no long range correlations in the wavefunction the 
coefficient of the n 2 term is simply the square of the groundstate energy per site. 
Of more interest is the linear term, which gives the variance per site, 

((H - nE ) 2 ) n = (H 2 ) n - n 2 E = no 2 , (27) 

where Eq is the energy per site. The variance gives a simple and reliable mea- 
sure of how close a given wavefunction is to being an eigenstate of H. In DMRG 
calculations the truncation error is usually used for this purpose Q], however the 
truncation error is not an intrinsic property of the wavefunction itself but rather is 
a byproduct of the particular choice of algorithm. Indeed, for some variants of MPS 
algorithms, such as single-site DMRG and TEBD in the limit of small time-step, 
the truncation error is identically zero and is therefore not useful. The variance 
however is an observable that can be used on any MPS wavefunction, irrespective 
of how the state was originally obtained. To test the utility of this approach, we 
have calculated the variance per site in the infinite size limit of the isotropic spin 
1/2 Heisenberg model, 

H = ^] (—SfSj - SfSj + SfSjj , (28) 

<ij> 

using iDMRG with a 1-site unit cell. As the number of states is varied, the variance 
and energy change according to Fig. [TJ Similarly to the well-known case of the trun- 
cation error in DMRG [J] , the variational energy is a linear function of the variance 
a 2 . For this example a linear fit using the 10 most accurate data points gives a 
groundstate energy of —0.443147181, correct to 9 significant figures compared with 
the exact result Q of 1/4 - In 2 ~ -0.4431471805599 .... By comparison, the best 
variational result in the calculated data, for 200 states kept, is —0.44314711, correct 
to only 7 significant figures. The accuracy of the fit is well captured by the standard 
error of the fit, of ue = 2 x 10 -9 , whereas linear fits calculated via the truncation 
error give an error that is often somewhat too small. This can be explained as a 
systematic error in the DMRG algorithm, by the convergence as the number of ba- 
sis states is changed where the truncation error typically changes quite smoothly. 
The truncation error measures convergence of the eigenvalues of the reduced den- 
sity matrix, which can be misleading because the eigenvalues converge much faster 
than the eigenvectors. This is the critical difference to the variance, which instead 
measures the convergence of the eigenvectors themselves and is therefore a more 
robust measure of convergence. 

5. Summary and conclusions 

In this paper, we have developed an algorithm for obtaining the expectation value 
and corresponding block operators of arbitrary triangular (Schur) matrix product 
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operators, which are obtained as a polynomial function of the lattice size n, in the 
large n limit. Operators of this type are useful for many purposes, and we have 
described a simple example where the square of the Hamiltonian can be used as 
an effective convergence test in numerical calculations, superseding the well-known 
truncation error from DMRG. In [20], one of us showed that for finite size MPS, it 
is practical to calculate power series and perturbative expansions to a dozen or so 
orders, and the algorithm we have presented here is similarly useful for obtaining 
the first few terms of an expansion, for example to obtain the pole expansion of the 
Green's function G(w,k). Other applications include obtaining expectation values 
of operators at finite momentum directly in the thermodynamic limit, which share 
a similar structure to the ansatz for an excited state studied previously by Ostlund 
and Rommer [0]. We expect that this algorithm will be an important component 
for constructing improved variational algorithms based around similar excited state 
ansatze. 
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